/*	This program produces Figure 3 */

***** Set directories 
local dir_clean 	"~/Dropbox/Retirement gaming/clean"
local dir_output 	"~/Dropbox/Retirement gaming/output/dataverse"


use "`dir_clean'/ech_clean.dta", clear

global listcontrols = " i.t i.fsize i.ciiu2  i.year#i.fsize i.year#i.ciiu2   " 


** SELF EMPLOYED EARNINGS (Figure 3a - earnings)
preserve
keep if self_empl==1
*Normalize hours
sum W if age_centered==-1
local meanw=r(mean)
g Wz=W/`meanw'
* Age dummies
tab age_centered, gen(age_centered_dums)
drop age_centered_dums5 
* Regression
reghdfe Wz age_centered_dums* [aw=weight], absorb(${listcontrols}) vce(robust)
* Sample size
global n=e(N)
* Coefficients and sd
for any beta sd: gen X = .
forvalues d = 1(1)13 {
	capture qui replace beta = _b[age_centered_dums`d'] if age_centered==`d'-6
	capture qui replace sd   = _se[age_centered_dums`d'] if age_centered==`d'-6
}
replace beta=0 if age_centered==-1
replace sd=0 if age_centered==-1
* Plot the results
collapse beta* sd*, by(age_centered)
gen sd_top = beta + 1.96*sd	
gen sd_bot = beta - 1.96*sd
gen sd_top10 = beta + 1.645*sd	
gen sd_bot10 = beta - 1.645*sd
sum sd*
twoway rcap sd_top sd_bot age_centered, ///
	cmissing(n) lwidth(thin) lpattern(solid) xsc(r(-5 7)) xlabel(-5(1)7)  xline(-1) xline(3.5, lp(dash) ) ysc(r(-0.15 .15)) ylabel(-0.15(0.05).15,grid) ///
	|| rcap sd_top10 sd_bot10 age_centered, cmissing(n) lwidth(thin)  lpattern(solid) xsc(r(-5 7)) ///
	|| scatter beta age_centered, scheme(s1mono) lpattern(solid) yline(0, lcolor(gs2)) ///
	title("Earnings, relative to age-49 average")  note("N=${n}") ///
	legend(off) ytitle("Estimated Coefficients") xtitle("Years Relative to Age 50")   
graph export "`dir_output'/figure3a_earn.png",  replace 	
restore


** SELF EMPLOYED HOURS (Figure 3a - hours)
preserve
keep if self_empl==1
*Normalize hours
sum hours if age_centered==-1
local meanh=r(mean)
g hoursz=hours/`meanh'
* Age dummies
tab age_centered, gen(age_centered_dums)
drop age_centered_dums5 // drop year before ref age, so everything becomes relative to that
* Regression
reghdfe hoursz age_centered_dums* [aw=weight], absorb(${listcontrols}) vce(robust)
* Sample size
global n=e(N)
* Coefficients and sd
for any beta sd: gen X = .
forvalues d = 1(1)13 {
	capture qui replace beta = _b[age_centered_dums`d'] if age_centered==`d'-6
	capture qui replace sd   = _se[age_centered_dums`d'] if age_centered==`d'-6
}
replace beta=0 if age_centered==-1
replace sd=0 if age_centered==-1
* Plot the results
collapse beta* sd*, by(age_centered)
gen sd_top = beta + 1.96*sd	
gen sd_bot = beta - 1.96*sd
gen sd_top10 = beta + 1.645*sd	
gen sd_bot10 = beta - 1.645*sd
sum sd*
twoway rcap sd_top sd_bot age_centered, ///
	cmissing(n) lwidth(thin)  lpattern(solid) xsc(r(-5 7)) xlabel(-5(1)7)  xline(-1) xline(3.5, lp(dash) ) ysc(r(-0.1 .1)) ylabel(-0.1(0.05).1,grid)   ///
	|| rcap sd_top10 sd_bot10 age_centered, cmissing(n) lwidth(thin) lpattern(solid) xsc(r(-5 7)) ///
	|| scatter beta age_centered,  scheme(s1mono)  lpattern(solid) yline(0, lcolor(gs2)) ///
	title("Hours of work, relative to age-49 average")  note("N=${n}") ///
	legend(off) ytitle("Estimated Coefficients") xtitle("Years Relative to Age 50")   
graph export "`dir_output'/figure3a_hours.png",  replace 	
restore

	

** EMPLOYED FIRM<10 EARNINGS
preserve
keep if empl==1 & small==1
*Normalize wages
sum W if age_centered==-1
local meanw=r(mean)
g Wz=W/`meanw'
* Age dummies
tab age_centered, gen(age_centered_dums)
drop age_centered_dums5 
* Regression
reghdfe Wz age_centered_dums*  [aw=weight], absorb(${listcontrols}) vce(robust)
* Sample size
global n=e(N)
* Coefficients and sd
for any beta sd: gen X = .
forvalues d = 1(1)13 {
	capture qui replace beta = _b[age_centered_dums`d'] if age_centered==`d'-6
	capture qui replace sd   = _se[age_centered_dums`d'] if age_centered==`d'-6
}
replace beta=0 if age_centered==-1
replace sd=0 if age_centered==-1
* Plot the results
collapse beta* sd*, by(age_centered)
gen sd_top = beta + 1.96*sd	
gen sd_bot = beta - 1.96*sd
gen sd_top10 = beta + 1.645*sd	
gen sd_bot10 = beta - 1.645*sd
sum sd*
twoway rcap sd_top sd_bot age_centered, ///
	cmissing(n) lwidth(thin)  lpattern(solid) xsc(r(-5 7)) xlabel(-5(1)7)  xline(-1) xline(3.5, lp(dash) ) ysc(r(-0.15 .15)) ylabel(-0.15(0.05).15,grid)  ///
	|| rcap sd_top10 sd_bot10 age_centered, cmissing(n) lwidth(thin) lpattern(solid) xsc(r(-5 7)) ///
	|| scatter beta age_centered,  scheme(s1mono) lpattern(solid) yline(0, lcolor(gs2)) ///
	title("Earnings, relative to age-49 average") note("N=${n}") ///
	legend(off) ytitle("Estimated Coefficients") xtitle("Years Relative to Age 50")   
graph export "`dir_output'/figure3b_earn.png",  replace 	
restore


** EMPLOYED FIRM<10 HOURS
preserve
keep if empl==1 & small==1
*Normalize hours
sum hours if age_centered==-1
local meanh=r(mean)
g hoursz=hours/`meanh'		
* Age dummies
tab age_centered, gen(age_centered_dums)
drop age_centered_dums5 
* Regression
reghdfe hoursz age_centered_dums* [aw=weight], absorb(${listcontrols}) vce(robust)
* Sample size
global n=e(N)
* Coefficients and sd
for any beta sd: gen X = .
forvalues d = 1(1)13 {
	capture qui replace beta = _b[age_centered_dums`d'] if age_centered==`d'-6
	capture qui replace sd   = _se[age_centered_dums`d'] if age_centered==`d'-6
}
replace beta=0 if age_centered==-1
replace sd=0 if age_centered==-1
* Plot the results
collapse beta* sd*, by(age_centered)
gen sd_top = beta + 1.96*sd	
gen sd_bot = beta - 1.96*sd
gen sd_top10 = beta + 1.645*sd	
gen sd_bot10 = beta - 1.645*sd
sum sd*
twoway rcap sd_top sd_bot age_centered, ///
	cmissing(n) lwidth(thin)  lpattern(solid) xsc(r(-5 7)) xlabel(-5(1)7)  xline(-1) xline(3.5, lp(dash) ) ysc(r(-0.1 .1)) ylabel(-0.1(0.05).1,grid) ///
	|| rcap sd_top10 sd_bot10 age_centered, cmissing(n) lwidth(thin) lpattern(solid) xsc(r(-5 7)) ///
	|| scatter beta age_centered, scheme(s1mono)  lpattern(solid) yline(0, lcolor(gs2)) ///
	title("Hours of work, relative to age-49 average")  note("N=${n}") ///
	legend(off) ytitle("Estimated Coefficients") xtitle("Years Relative to Age 50")   
graph export "`dir_output'/figure3b_hours.png",  replace 	
restore




** EMPLOYED FIRM>=10 EARNINGS
preserve
keep if empl==1 & small==0 
*Normalize wages
sum W if age_centered==-1
local meanw=r(mean)
g Wz=W/`meanw'
* Age dummies
tab age_centered, gen(age_centered_dums)
drop age_centered_dums5 
* Regression
reghdfe Wz age_centered_dums* [aw=weight], absorb(${listcontrols}) vce(robust)
* Sample size
global n=e(N)
* Coefficients and sd
for any beta sd: gen X = .
forvalues d = 1(1)13 {
	capture qui replace beta = _b[age_centered_dums`d'] if age_centered==`d'-6
	capture qui replace sd   = _se[age_centered_dums`d'] if age_centered==`d'-6
}
replace beta=0 if age_centered==-1
replace sd=0 if age_centered==-1
* Plot the results
collapse beta* sd*, by(age_centered)
gen sd_top = beta + 1.96*sd	
gen sd_bot = beta - 1.96*sd
gen sd_top10 = beta + 1.645*sd	
gen sd_bot10 = beta - 1.645*sd
sum sd*
twoway rcap sd_top sd_bot age_centered, ///
	cmissing(n) lwidth(thin)  lpattern(solid) xsc(r(-5 7)) xlabel(-5(1)7)  xline(-1) xline(3.5, lp(dash) ) ysc(r(-0.15 .15)) ylabel(-0.15(0.05).15,grid)   ///
	|| rcap sd_top10 sd_bot10 age_centered, cmissing(n) lwidth(thin) lpattern(solid) xsc(r(-5 7)) ///
	|| scatter beta age_centered, scheme(s1mono)  lpattern(solid) yline(0, lcolor(gs2)) ///
	title("Earnings, relative to age-49 average")  note("N=${n}") ///
	legend(off) ytitle("Estimated Coefficients") xtitle("Years Relative to Age 50")   
graph export "`dir_output'/figure3c_earn.png",  replace 	
restore


** EMPLOYED FIRM>=10 HOURS
preserve
keep if empl==1 & small==0 
*Normalize hours
sum hours if age_centered==-1
local meanh=r(mean)
g hoursz=hours/`meanh'
* Age dummies
tab age_centered, gen(age_centered_dums)
drop age_centered_dums5 
* Regression
reghdfe hoursz age_centered_dums* [aw=weight], absorb(${listcontrols}) vce(robust)
* Sample size
global n=e(N)
* Coefficients and sd
for any beta sd: gen X = .
forvalues d = 1(1)13 {
	capture qui replace beta = _b[age_centered_dums`d'] if age_centered==`d'-6
	capture qui replace sd   = _se[age_centered_dums`d'] if age_centered==`d'-6
}
replace beta=0 if age_centered==-1
replace sd=0 if age_centered==-1
* Plot the results
collapse beta* sd*, by(age_centered)
gen sd_top = beta + 1.96*sd	
gen sd_bot = beta - 1.96*sd
gen sd_top10 = beta + 1.645*sd	
gen sd_bot10 = beta - 1.645*sd
sum sd*
twoway rcap sd_top sd_bot age_centered, ///
	cmissing(n) lwidth(thin)  lpattern(solid) xsc(r(-5 7)) xlabel(-5(1)7)  xline(-1) xline(3.5, lp(dash) ) ysc(r(-0.1 .1)) ylabel(-0.1(0.05).1,grid) ///
	|| rcap sd_top10 sd_bot10 age_centered, cmissing(n) lwidth(thin) lpattern(solid) xsc(r(-5 7)) ///
	|| scatter beta age_centered, scheme(s1mono)  lpattern(solid) yline(0, lcolor(gs2)) ///
	title("Hours of work, relative to age-49 average") note("N=${n}") ///
	legend(off) ytitle("Estimated Coefficients") xtitle("Years Relative to Age 50")   
graph export "`dir_output'/figure3c_hours.png",  replace 	
restore





clear all
exit
 
